MODULE GLOBVAR
  USE NRUTIL, ONLY : ARRAY1D, IARRAY1D, ARRAY2D, IARRAY2D, ARRAY3D, IARRAY3D, ARRAY7D, IARRAY7D, NR_BIG
  IMPLICIT NONE

  real(8), parameter :: &
        & s_tol = 1.0d-10, &
        & Euler_Constant = 0.5772d0, &
        & gr_exitcode = -1.0d6

  integer, parameter :: &
        & Icol = 0,  &  !High School (Icol==0) or College (Icol==1)
        & IHasHealthAge = min(0, 50),   &  !0-no health, 50-has different health status starting from age 50
        & Ipt = 0,   &       !Ipt=0: Full time only, Ipt=1: Part time and Full time. 
        & Ipt_coef = 1, &    ! 1-polynomial, centered at N_0; 2-linear + alpha; 3-centered at 45;
                             ! 15--N_plpt=4, search center; 21: 1+FS, 22: 2+FS
        & N_plpt = 3 + Ipt_coef/10,  &     !3 dimension of plpt_V
        & N_phi_t = 3, &  !consumption shifter: t, t^2, or t^3
        !& N_phi_t = 0, &  !cs_noage*: consumption shifter has no age curve: exp(0) 
        & N_ssa_u = 3, &     !3: pssa_62, pssa_65, pssa_65_t
        !& N_ssa_u = 0, &    !ssa_noage*: pssa_62, pssa_65, pssa_65_t are all 0.0d0                
        & N_pbequest = 1, &  !estimate pb1 only and set pb2 = 222.0d0 !as in French and Jones (2011, Econometrica)
        !& N_pbequest = 2, & !estimate both pb1 and pb2
        & N_depreciation = 1,  &  ! 1--H depreciates at the same rate on/off job
        !& N_depreciation = 2,  &  ! 1--H depreciates at different rate on or off job        
        & IC_L_nonseparable = 0,  &  !0-separable, 1-nonseparable
        & gcI_IncreaseLE = 0, &  !increase life expectancy by this much
        & N_0 = 18 + Icol*4,    &         !period starts
        & N_T = 80 + gcI_IncreaseLE,  &   !max age. Life expectancy for white male in US is 74.13 
                                          !in 2000 and 76.5 in 2010.
        & N_T_grid = N_T - N_0 + 1,   &   !t: 1:N_T_grid in state grids.
        & N_FS = 3,  &  !family status: 1 (single or divoiced), 2 (married spouse not working), 3 (married spouse working)
        & N_FS_spinc_int = 2, &  !points for Gauss-Hermit integration over spouse income
        & N_data_0 = 22 + Icol*4 ,   &   !where data begins
        & N_data_10 = 65,  &   !where data ends. For SIPP
        & N_data_0_health = 50, &  !where lfpr diff in health starts, from 50
        & N_data_0_c = N_data_0,  &  !
        & N_data_0_ss = 62, N_data_10_ss = 70,  &
        & N_data_0_pt = N_0, N_data_10_pt = 65,  &  !for part time
        & N_moments_type = 6+min(3,IHasHealthAge) + Ipt,  &     ! number of types of all moments 
                                     ! 1 lfpr, 2 lnw, 3 lnw_sd, 4 lnw_fe, 5 c, 6 ss, 
                                     ! (if with Health) 7 diff_lfpr_E2G, 8 diff_lfpr_G2B, 9 diff_lfpr_B2D, 7/10 lfpr_pt
        & N_moments_all = (N_data_10-N_data_0+1)*4    &  !1 lfpr, 2 lnw, 3 lnw_sd, 4 lnw_fe
                   &  + (N_data_10-N_data_0_c+1)      &  !5 consumption
                   &  + (N_data_10_ss-N_data_0_ss+1)  &  !6 ss
                   &  + min(3,IHasHealthAge)*(N_data_10-N_data_0_health+1)  &  !7 diff_lfpr_e2g, 8 diff_lfpr_g2b, 9 diff_lfpr_g2b
                   &  + Ipt * (N_data_10_pt-N_data_0_pt+1)  &   !(7 or 10 Part Time)
                   &  + 2 + 1,  &  !lfpr transition (2), and wage depreciation after one period of unemployment (1)
        & N_XI_int = 5, &      !points for Gauss-Hermit integration over XI
        & N_het_a0 = 3, N_het_pi = 3, N_het_H0 = 2, &   ! Number of heterogeneity dimensions. N_het_H0 takes 1 or more.
        & N_PAR = 13+N_pbequest+N_phi_t+5*min(1,IHasHealthAge)+(min(N_het_a0,2)-1)*(min(N_het_pi,2)-1)*(3+2*(min(N_het_H0,2)-1)) &
        &           +min(1,IC_L_nonseparable) + N_ssa_u +Ipt*N_plpt + (N_depreciation-1), &
        & N_PAR_mpi = N_PAR + 2, &   !N_PAR+1/2 are for controls
        & N_Sim = 20000 + 30000 * min(1, (Ipt+min(1,IHasHealthAge))) - 10000 * min(1, (Ipt*IHasHealthAge))
        !50k for Part Time PT or Health; adjust to 40K for PTHealth
        ! 40K is sufficient for pthealth. 80K generates very similar results.
  
  logical :: gc_nogridsearch, gcL_debugprint, gcL_noFS, gcL_partrans,  &
          &  gcL_spinc_discrete, gcL_XI_discrete, &
          &  gcL_printvalue, gcL_cal_shadowprice
  integer(4) :: N_A, VN_A(N_0:N_T), &   !asset space
        & N_H, &    !human capital space
        & N_AIME, VN_AIME(N_0:N_T), &    !AIME space
        & N_HEALTH, &  !Health
        & N_RECSS,  &    !ret or not, taking ss or not
        & N_Data         !Number of Data observations
        
  !===========================================================
  !parameters, starting with p
  !===========================================================
  real(8) :: pa0, pa_health1E, pa_health3B, pa_health3B_t, pa_health4D, pa_health4D_t, plpt_V(N_plpt), plpt_t(N_0:N_T,N_FS),  & 
          &  pa_epsilon, palpha_I, palpha_H, psigma, psigma_onjob, psigma_offjob, ppi,  &
          &  pssa_62, pssa_65, pssa_65_t,    &
          &  pXI_sd, pinitHmean, pinitHsd, pinitH_a0, pinitH_pi, &  !parameters to be estimated, including initial distributions.
          &  pphi(N_phi_t+2), petac, pb1, pa1, pa2, pinitA, grinitAsd
  real(8) :: pr, pr_inv, pw(100), pbeta, pb2, gr_XI_sd,gr_XI_mean  !consumption, bequest, H shock
          ! parameter vector and domain vector
  real(8) :: gr_paradomain(N_PAR,2), gr_Vpar(N_PAR), gr_Vpar_mpi(N_PAR_mpi), gr_Vpar_orig(N_PAR),  &
           & p1_Vhet_x(N_het_a0,2), p2_Vhet_x(N_het_pi,2),  &  !discretized parameters of (pa0,pi)
           & pinitA_domain(1,2), pintA_Vmpi(1:2)  !domain of 
  integer(4) :: gI_het_group_p1(N_het_a0 * N_het_pi), gI_het_group_p2(N_het_a0 * N_het_pi), &  !mapping from group to index of (pa0,pi)
          &  gI_het_pos(N_het_a0, N_het_pi)  !mapping from index of (pa0,pi) to group
 
  real(8) :: gr_bequest_AP(N_FS,0:(Ipt+1)),  &
           & gr_LBD_t, gr_LBD_tsq, &
           & grTotHours, grCF,  &  !total hours, for normalization; consumption floor
           & grUVmultiplier, grUVmultiplier_ln, &
           & gr_sstax_1,gr_sstax_2,gr_sstax_3,gr_sse_cap, &
           & gr_psi(N_0:N_T,N_FS,0:(Ipt+1)), &
           & gr_atilde(3),    &
           & gr_leisure_taste(N_0:N_T,1:4,N_FS),    &  !(git,giHEALTH,giFS)
           & gr_leisure_condexp(N_0:N_T,1:4,N_FS),    &
           & grAdultEquivScale(5)
    
  real(8) :: gr_XI_xw(N_XI_int,3),  &  !1-x (standard normal), 2-weights, 3-actual XI
          &  gr_FS_spinc_xw(N_0:N_T, N_FS_spinc_int,3),  &
          &  gr_health_matrix(N_0:N_T,1:4,0:4)  !1 age, 2 status (1exce 2good 3bad 4DI), 3 (0dist 1to exce 2to good 3to bad 4toDI) !for health transition
 
  !for family transition and spouse income
  integer(4) :: giFS   !family or marital status
  real(8) :: gr_FS_trans(N_0:N_T,1:N_FS,1:N_FS)   !t, current, next
  real(8) :: gr_FS_inc(N_0:N_T,2)  !t, (mean, sd)
  
  !===========================================================
  !  policy related global variables
  !===========================================================
  integer(4) :: pERA, pNRA      ! early retirement age (no ret/ss before this age), normal retirement age
  real(8) :: pPIAbp(2), pETESTbp(2)  !PIA bend points, Earnings test bending point
  real(8) :: pPIAratio(3)            !PIA raito in Equation C4: 0.9, 0.32, 0.15 in the data
  real(8) :: gr_ftax_pbk(N_0:N_T, 8), gr_ftax_MPTR(N_0:N_T, 8), gr_ftax_pbase(N_0:N_T, 8)  !federal tax
  real(8) :: gr_SSDI_SGA, gr_SSI  !SSI benefit
  
  !state variables
  integer(4) :: git, git_pos, giA, giH, giAIME, giHEALTH, giRECSS, giRECSS_next    !indexing
  real(8) :: gA, gH, gAIME
  
  ! Grid
  TYPE(ARRAY1D) :: sgridcom, gr_Amax, sgridAnrlb, sgridAnrlb_BC  !state space grid, Asset borrowing lower bound--same as sgridA%(:,1)
  TYPE(ARRAY2D) :: sgridH, sgridAIME, sgridA, sgridA_BC
  TYPE(IARRAY1D) :: sgridRECSS  ! RECSS: 0-no ret / no ss; 1-ret/ss.  HEALTH: 1-good, 2-bad
  integer(4) :: gI_HgridIncr, gI_HgridFullNodes
      
  ! value function and policy functions
  ! (0) is for the case with leisure=0, (1) is for leisure=1 (not work)
  ! _EV is the Emax
  TYPE(ARRAY7D) :: dvV(0:(Ipt+1)), dvV_EV, &  !value function: A, AIME, HEALTH, RECSS, H, t
                 & dcc(0:(Ipt+1)),   &        !policy (control) functions: consumption,
                 & dcI(0:Ipt), depsilon(0:Ipt)        !investment, threshold of working or not; 0 for Leisure=0 and 1 for Leisure=pt
  TYPE(IARRAY7D) :: dciRECSS(0:(Ipt+1))       !policy function of ret/ss
  
  !for Investment-Work decomposition
  logical :: gcL_CF_IW_decompose
  integer :: gcI_CF_IW
  TYPE(ARRAY7D) :: dvV0(0:(Ipt+1)), dvV_EV0, &  !value function: A, AIME, HEALTH, RECSS, H, t
                 & dcc0(0:(Ipt+1)),   &        !policy (control) functions: consumption,
                 & dcI0(0:Ipt), depsilon0(0:Ipt)        !investment, threshold of working or not; 0 for Leisure=0 and 1 for Leisure=pt
    
  !global controls
  integer(4) :: gcI_itercount
  logical :: gcL_Iscapped, gcL_IsDRC, gcL_IsOutput,  &
          &  gcL_cI_Simplex, gcL_triangle, gcL_searchinitA,  &
          &  gcL_searchptonly, gcL_skippt,  &   !skip pt when calculating lfpr, as in data
          &  gcL_searchSSAonly,  &   !search SSA parameters only
          &  gcL_searchBC   !search for the ratio of borrowing constraint
          
  integer(4) :: gcI_robust, gcI_elas_t, gcI_taxhike_t, gcI_CF, gI_p_L
  integer(4) :: gcI_imax, gcI_inte_c_log, gcI_inte_v_log, gcI_inte_IL_log !0-normal, 1-x and y log
  real(8) :: gcF_vtol, grMinlnw, grMinDistance
  !gcF_vtol is for interior solutions
  !dp_solve use a different one which is set 1.0d-10 in the code, but it is still very fast.
  
  !used in Module cmaxee
  real(8) :: gr_p_Ibound(0:2,2), gr_p_cbound(2), gr_p_I, gr_p_L, gr_cc, gr_BL, gr_LowerBound, gr_BC_Ratio

  !=================== data and simulation ==========================================
  TYPE(ARRAY2D) :: MSM_Moments_Data !data, age. 
  TYPE(ARRAY1D) :: cps_aimeshare !share of min wage at AIME
  real(8) :: MSM_Moments_Data_ss(N_data_0_ss:N_data_10_ss,2)    !1 age, 2 (1)ssa (2)irecss
  
  ! output
  TYPE(ARRAY3D) :: MSM_Indiv, MSM_Indiv_Alt 
                  !First dimension: 1-A,2-H,3-AIME,4-c,5-I,6-L,7-Labor,8-Value,9-leisure error
                  !Second dimension: from N_0 to N_T; Third dimension: 1-N_SIM
  TYPE(ARRAY2D) :: MSM_Indiv_initX !Initial condition at N_0 --- {N_SIM, 1A 2H} 
  TYPE(IARRAY2D) :: MSM_Indiv_Ihet  !N_Sim; second dimension: 1. heterogeneity type, 2-observed first year, 3-observed last year
  TYPE(IARRAY2D) :: MSM_Indiv_Irecss, & ! receiving ss (1) or not (0). First dimension: 1-N_SIM; Second dimension: from N_0 to N_T. 
                  & MSM_Indiv_Ihealth, MSM_Indiv_IFS
  TYPE(IARRAY3D) :: MSM_Indiv_Irecss_Alt ! receiving ss (1) or not (0).
  
  TYPE(ARRAY2D) :: MSM_Moments_Sim !Simulated Method of Moments
  TYPE(ARRAY2D) :: MSM_wgt_scale, gmmse_S, MSM_wgt !rescale moments, weight in the GMM at MSM
    
  TYPE(IARRAY2D) :: MSM_sim4m_i !index of hrs_sim4m
  !==============================================================================================================================
    
  
  !===================================================================
  ! Variables for MPI
  integer(4) :: myid_mpi, numprocs_mpi, master_mpi,  &
          &  gI_vsize_mpi, gI_Indiv_size_mpi, gI_Indiv_Alt_size_mpi,  &
          &  new_comm_mpi, new_id_mpi, new_numprocs_mpi, color_mpi, gI_displacement_mpi
  TYPE(IARRAY2D) :: gI_VTasks_mpi, gI_Indiv_Tasks_mpi, gI_Vcolor
  !===================================================================
  
  CHARACTER(len=255) :: cwd, cwd_indiv
  character(1) :: gc_filenamehet
  integer(4) :: gI_DT0(8), gI_DT1(8)
  
  
END MODULE GLOBVAR
